' todo:
' megaparticles by nvidia
' convert to cuda
' better water code sph. or something. this takes lots of cpu :(
' pad grid so corner pixels show.... or just clip window?


' to interpolate grid or not. smoother surfs at lower res.

SuperStrict

Framework BRL.GLMax2D
Import BRL.PngLoader
Import BRL.Standardio
Import Pub.GLEW
Import BRL.Random
SetGraphicsDriver GLMax2DDriver()
Global xres#=800
Global yres#=600

'-----Grid
Global grid_xres:Int = 80*3
Global grid_yres:Int = 60*3
Global cell_w# = xres/Float(grid_xres)
Global cell_h# = yres/Float(grid_yres)

Global grid:TBlock[grid_xres+2,grid_yres+2] 'todo: +2 units of padding in for marching squares bounds.

For Local i:Int=0 Until grid.dimensions()[0]
For Local j:Int=0 Until grid.dimensions()[1]

	grid[i,j] = New TBlock	
Next
Next

Global msquares:TMarchingSquares = New TMarchingSquares  
'------

Graphics xres,yres,0,0,2


'SetMaskColor(255,255,255)
SetMaskColor(0,0,0)
Global texParticle:Int = GLTexFromPixmap( LoadPixmapPNG("metaball.png") )'GLTexFromPixmap

Global water:waterparticle[] = New waterparticle[3000]
For Local i:Int = 0 Until Len water
	water[i] = New waterparticle
Next

Global num_particles:Int = 0
Global radi# = 10

Global isoMinimum:Int = 0

Global doAnimation:Int = False
Global doMarcingSquares:Int = False

initPointSprite()


Global e# = 0.9 'restitution


Local tim:Int, oldtim:Int
Local fps:Int = 0
Local fpsTimer:Int = 0
Local numframes:Int = 0
Global dt# = 0

While Not KeyHit(KEY_ESCAPE) And Not AppTerminate()
	tim = MilliSecs()
	dt = (tim-oldtim) / 1000.0
	oldtim = tim
	
	If tim > fpsTimer
		fps = numframes
		numframes = 0
		fpsTimer = tim + 1000
	EndIf
	numframes:+1
	
	'clear grid
	For Local b:TBlock = EachIn grid
		b.num_particles = 0
	Next

	Local perfTim1:Int = MilliSecs()
	dynamics()
	Local perfTim2:Int = MilliSecs()
	collisions() 'cant pause collisions,or else we get a bad grid pos
	Local perfTim3:Int = MilliSecs()

	Local first:waterparticle = water[0]
	If MouseDown(2)
	If first <> Null
		'todo: find a way to give momentum
		'first.vx = MouseX() - first.px
		'first.vy = MouseY() - first.py
		first.px = MouseX()
		first.py = MouseY()
	EndIf
	EndIf
	
	If MouseDown(3) And num_particles > 0 num_particles:-1
	
	If MouseDown(1) And num_particles < Len(water)
		Local activatedParticle:waterparticle = water[num_particles]
		activatedParticle.px = MouseX() + Rand(-1,1)
		activatedParticle.py = MouseY() + Rand(-1,1)
		activatedParticle.nx = activatedParticle.px
		activatedParticle.ny = activatedParticle.py
		activatedParticle.ay = 0
		num_particles:+1
	EndIf
	
	If KeyHit(KEY_SPACE)
		doAnimation= Not doAnimation
	EndIf

	If KeyHit(KEY_R)
		num_particles = 0
	EndIf

	If KeyHit(KEY_M)
		doMarcingSquares = Not doMarcingSquares
	EndIf
	

	If doMarcingSquares
		glPointSize(8)
		
		glEnable(GL_TEXTURE_2D)
		glBindTexture(GL_TEXTURE_2D, texParticle)
		glEnable(GL_BLEND)
		glBlendFunc(GL_ONE, GL_ONE) 'Additive blending
		glColor4f(0.1,0.2,0.8,.5)
		glEnable(GL_POINT_SPRITE)
		glBegin(GL_POINTS)
			msquares.getOutline(grid)
		glEnd()
		glDisable(GL_POINT_SPRITE)	
		
	Else
		renderworld()
	EndIf
	Local perfTim4:Int = MilliSecs()
	
	glPointSize(1)
	'drawgrid()
	
	SetAlpha 1
	SetColor 255,255,255
	SetBlend MASKBLEND
	DrawText "dt "+dt, 0, 15
	DrawText "# "+num_particles, 0, 15*2
	DrawText "fps: "+fps, 0, 15*3
	DrawText "keys: R, M, space, mouse 1,2,3", 0, 15*4
	
	DrawText "dynamics : "+(perftim2-perftim1), 0, 15*6
	DrawText "collision: "+(perftim3-perftim2), 0, 15*7
	DrawText "render   : "+(perftim4-perftim3), 0, 15*8

	Flip
	Cls
	
	?debug
	checkGLError()
	?
Wend

Function renderworld()
	glEnable(GL_TEXTURE_2D)
	glBindTexture(GL_TEXTURE_2D, texParticle)

	glEnable(GL_BLEND)
	glBlendFunc(GL_ONE, GL_ONE) 'Additive blending
	'glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA ) 'additive without overblending

	glPointSize radi*4'*4.0
	
	glColor4f(0.1,0.2,0.8,.6)

	glEnable(GL_POINT_SPRITE)
	glBegin(GL_POINTS)
		
	Local i:Int
	For i = 0 Until num_particles
		Local w:waterparticle = water[i]
		glVertex2f(w.px,w.py)
	Next
	
	glEnd()
	glDisable(GL_POINT_SPRITE)
End Function

Function collisions()
	Local i:Int
	Local j:Int
	
	For i=0 Until num_particles
		Local a:waterparticle = water[i]
		For j=i+1 Until num_particles
			Local b:waterparticle = water[j]
			
			Local dx# = (b.px - a.px)
			Local dy# = (b.py - a.py)
			Local distSq# = dx*dx+dy*dy
			
			Local vxa# = 0.0; Local vya# = 0.0
			Local vxb# = 0.0; Local vyb# = 0.0
			
			If distSq < radi*radi*4 'all balls same size. 2 balls intersect
				Local dist# = Sqr(distSq)
				vxa = a.vx; vya = a.vy;
				vxb = b.vx; vyb = b.vy;
				Local nx# = b.px - a.px
				Local ny# = b.py - a.py
				
				If Abs(dist) > 1e-5
					Local rdist# = 1/dist
					Local cdx# = nx * rdist
					Local cdy# = ny * rdist
					
					Local dd# = 2*radi-dist
					'Move apart by intersection dist
					a.nx :- dd*cdx
					a.ny :- dd*cdy
					b.nx :+ dd*cdx
					b.ny :+ dd*cdy
					
					
					'New velocities
					'Local am# = Sqr( a.vx*a.vx + a.vy*a.vy )
					Local am# = ( a.vx*a.vx + a.vy*a.vy )


					If Abs(am) > 1e-5*1e-5
						'dot prod mult by restitution
						Local cosam# = cdx * a.vx + cdy * a.vy * e
						vxa :- cdx*cosam
						vya :- cdy*cosam
						vxb :+ cdx*cosam
						vyb :+ cdy*cosam
					EndIf
					
					'Local bm# = Sqr( b.vx*b.vx + b.vy*b.vy )
					Local bm# = ( b.vx*b.vx + b.vy*b.vy )


					If Abs(bm) > 1e-5*1e-5
						'dot prod mult by restitution
						Local cosam# = cdx * b.vx + cdy * b.vy * e
						vxa :+ cdx*cosam
						vya :+ cdy*cosam
						vxb :- cdx*cosam
						vyb :- cdy*cosam
					EndIf
					
				EndIf
				
				a.vx = vxa
				a.vy = vya
				b.vx = vxb
				b.vy = vyb
			EndIf
				

		Next
	Next
End Function


Function dynamics()
	'Dynamics
	'chrystal-like behaviour. "stolen" from paul bourkes site.
	Local i:Int
	For i=0 Until num_particles
		Local w:waterparticle = water[i]

		
		w.px = w.nx
		w.py = w.ny
	
		Local tempx# = w.ax*dt
		w.px = w.px + dt*(w.vx+tempx*.5)
		w.vx = w.vx+ tempx

		Local tempy# = w.ay*dt
		w.py = w.py + dt*(w.vy+tempy*.5)
		w.vy = w.vy+ tempy
			
			
		If doAnimation 'do gravity really
			w.ax = 0
			w.ay = 100'200
		else
			w.ax = 0
			w.ay = 0'200	
		EndIf

		'Boundary conditions
		If w.px < w.r 
			w.px = w.r 
			w.vx = -w.vx * e
		EndIf
		
		If w.py < w.r 
			w.py = w.r 
			w.vy = -w.vy * e
		EndIf
		
		If w.px > xres-w.r
			w.px = xres-w.r 
			w.vx = -w.vx * e
		EndIf
		
		If w.py > yres-w.r 
			w.py = yres-w.r	 
			w.vy = -w.vy * e
		EndIf
		
		w.nx = w.px
		w.ny = w.py		
		
		'Assign to cell
		Local xcell:Int = Floor( w.px / xres * Float(grid_xres) )
		Local ycell:Int = Floor( w.py / yres * Float(grid_yres) )
		Assert xcell >= 0 And xcell < grid_xres
		Assert ycell >= 0 And ycell < grid_yres
		
		'grid[xcell,ycell].particles.addlast(w)
		grid[xcell,ycell].num_particles :+1
		'Up down, left right
		If xcell > 1 grid[xcell-1,ycell].num_particles :+1
		If xcell < grid_xres-1 grid[xcell+1,ycell].num_particles :+1
		If ycell > 1 grid[xcell,ycell-1].num_particles :+1
		If ycell < grid_yres-1 grid[xcell,ycell+1].num_particles :+1
		
		'Edges
		If xcell > 1 And xcell < grid_xres-1
		If ycell > 1 And ycell < grid_yres-1
			grid[xcell-1,ycell-1].num_particles :+1
			grid[xcell-1,ycell+1].num_particles :+1
			grid[xcell+1,ycell-1].num_particles :+1
			grid[xcell+1,ycell+1].num_particles :+1
		EndIf
		EndIf	
			
	Next
End Function

Function dist2D:Float(x1#,y1#,x2#,y2#)
	Return Sqr( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) )
End Function

Function distSquare:Float(x1#,y1#,x2#,y2#)
	Return (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)
End Function

Type waterparticle
	Field ox#, oy#
	Field px#, py#
	Field nx#, ny#
	Field vx#, vy#
	Field ax#, ay#
	Field r# = radi
End Type


Function checkGLError()
	Local checkedAllErrors:Int = False
	Local glerrorList$ = ""
	Local numErrors:Int = 0
	While checkedAllErrors=False
		Local errorCode:Int = glGetError()
		If errorCode = GL_NO_ERROR
			checkedAllErrors = True
		Else
			numErrors = numErrors + 1
			If errorCode = GL_INVALID_ENUM glerrorList:+" invalid enum "
			If errorCode = GL_INVALID_VALUE glerrorList:+" invalid value "
			If errorCode = GL_INVALID_OPERATION glerrorList:+" invalid operation "
			If errorCode = GL_STACK_OVERFLOW glerrorList:+" stack overflow "
			If errorCode = GL_STACK_UNDERFLOW glerrorList:+" stack underflow "
			If errorCode = GL_OUT_OF_MEMORY glerrorList:+" GL_OUT_OF_MEMORY "
			If errorCode = GL_TABLE_TOO_LARGE glerrorList:+" GL_TABLE_TOO_LARGE "
		EndIf
	Wend
	If glerrorList <> ""
		Print "#"+numErrors +" glGetErrors this frame: "+glerrorList
	EndIf
End Function

Function initPointSprite()
'Glew, opengl extension wrangler is used To see If an extension is supported.
	glewInit()	
	If(GL_POINT_SPRITE)
		Print "GL_POINT_SPRITE is supported"
		glEnable(GL_POINT_SPRITE)
	Else 
		Print "GL_POINT_SPRITE is Not supported"
		End
	EndIf
	
	glActiveTexture(GL_TEXTURE0);
	glEnable( GL_TEXTURE_2D );
	glTexEnvi(GL_POINT_SPRITE, GL_COORD_REPLACE, GL_TRUE); 'Generates tex coords for Point Sprite
	'modes: <a href="http://www.khronos.org/opengles/documentation/opengles1_0/html/glTexEnv.html" target="_blank">http://www.khronos.org/opengles/documentation/opengles1_0/html/glTexEnv.html</a>
	glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);  'GL_MODULATE, GL_DECAL, GL_BLEND, Or GL_REPLACE.
End Function

Type TBlock
	'Field particles:TList = New TList
	Field num_particles:Int = 0
End Type


Function drawrectopen(x0#,y0#,x1#,y1#)
	DrawLine(x0,y0,x1,y0)
	DrawLine(x1,y0,x1,y1)
	DrawLine(x1,y1,x0,y1)
	DrawLine(x0,y1,x0,y0)
End Function

Function drawGrid()
	SetColor(255,255,255)
	For Local i:Int=0 Until grid_xres
		For Local j:Int=0 Until grid_yres
			If grid[i,j].num_particles > 0 
				SetColor 0,0,255
			Else
				SetColor 255,255,255
			EndIf
			DrawRectOpen(i*cell_w, j*cell_h, (i+1)*cell_w, (j+1)*cell_h)
		Next
	Next
End Function 


Type TMarchingSquares
	Field bmd:TBlock[,]
	Field width:Int, height:Int
	Field variations_ccw:vec2[16]
	Field tracked:Int[,] 'used to remember if we've traced the pixel earlier
	Field points:TList 'list of traced points. Might contain dupes
	Field animated:Int = False
	
	Field foundLoop:Int = 0
	
	Method New()
		variations_ccw[ 0] = vec2.Create( 0, 0) ' illegal
		variations_ccw[ 1] = vec2.Create( 0, 1) ' d v
		variations_ccw[ 2] = vec2.Create(-1, 0) ' l <
		variations_ccw[ 3] = vec2.Create(-1, 0) ' l <
		variations_ccw[ 4] = vec2.Create( 1, 0) ' r >
		variations_ccw[ 5] = vec2.Create( 0, 1) ' d v
		variations_ccw[ 6] = vec2.Create(-1, 0) ' l <
		variations_ccw[ 7] = vec2.Create(-1, 0) ' l <
		variations_ccw[ 8] = vec2.Create( 0,-1) ' u ^
		variations_ccw[ 9] = vec2.Create( 0,-1) ' u ^
		variations_ccw[10] = vec2.Create( 0,-1) ' u ^
		variations_ccw[11] = vec2.Create( 0,-1) ' u ^
		variations_ccw[12] = vec2.Create( 1, 0) ' r >
		variations_ccw[13] = vec2.Create( 0, 1) ' d v
		variations_ccw[14] = vec2.Create( 1, 0) ' r >
		variations_ccw[15] = vec2.Create( 0, 0) ' illegal
	End Method
	
	
	Method getOutline:TList( bitmapdata:TBlock[,])
		bmd = bitmapdata
		width = grid_xres'bitmapdata.dimensions()[0]
		height = grid_yres'bitmapdata.dimensions()[1]
		
		Return searchEachPixel()
	End Method

	Method searchEachPixel:TList()
		tracked = New Int[width,height]
		points = New TList

		' reset error status
		foundLoop = 0
		
		'find pixels that have something in them and are untracked so far, run marching squares on pixel
		'keep boundries witin image (2x2) marching squares area 
		For Local y:Int=0 Until height
		For Local x:Int=0 Until width
			If(bmd[x+1,y+1].num_particles > isoMinimum) 'add 1 to x,y for padding
				If( Not tracked[x,y])
					marchingSquares(x,y);
				EndIf
			EndIf
		Next
		Next
	
		Return points
	End Method

  Method marchingSquares(x:Int,y:Int)
	Local pos:vec2 = vec2.Create( x-1, y-1 ) 'move back and up from pixel
	if pos.x < 0 or pos.y < 0 return
	Local start:vec2 = vec2.Create( pos.x, pos.y )
	Local outOfBounds:Int = pos.x < 0 Or pos.x > width-1 Or pos.y < 0 Or pos.y > height-1

	Assert Not outOfBounds
	
	Local iters:Int = 0
	
	While(True)
		iters:+1
		tracked[pos.x,pos.y] :+ 1
		Local dirTableIndex:Int = getNextEdgePoint(pos)
		outOfBounds = pos.x < 0 Or pos.x > width-1 Or pos.y < 0 Or pos.y > height-1

		If outOfBounds Then Return
		If( dirTableIndex=0 Or dirTableIndex=15 ) Then Return 'quit if illegal marching squares state
		If( pos.x=start.x And pos.y=start.y ) Then Return 'were back at start, finished trace
		If( tracked[pos.x,pos.y] > 2  )  'tracked this pixel too many times, prevent inf loops caused by bad pixel configurations.
			foundLoop:+1
			Return 
		EndIf
		
		
		If(1)
			glVertex2f(pos.x*cell_w, pos.y*cell_h)
			'points.addLast( vec2.Create(pos.x, pos.y) )
		Else
			Local px# = pos.x*cell_w
			Local py# = pos.y*cell_h
			DrawRectOpen(px,py, px+cell_w, py+cell_h)
		EndIf

		
		' Visualize marching squares progress
		If(animated)
			SetColor(0,50+5*iters Mod 150,0)
			Plot(pos.x,pos.y)
			
			If(iters Mod 5 = 0)  
				Flip()
				Delay(5)
			EndIf
		EndIf
		
	Wend
	
	End Method

	'use marching squares look-up-table to decide where to move next
	Method getNextEdgePoint:Int(pos:vec2)
		Local x:Int = pos.x+1; Local y:Int = pos.y+1 'add 1 to x,y for padding
		Local index:Int = 0
		If( bmd[x+1,y+1].num_particles > isoMinimum ) index :+ 1
		If( bmd[x  ,y+1].num_particles > isoMinimum ) index :+ 2
		If( bmd[x+1,y  ].num_particles > isoMinimum ) index :+ 4
		If( bmd[x  ,y  ].num_particles > isoMinimum ) index :+ 8	
		pos.x :+ variations_ccw[index].x; pos.y :+ variations_ccw[index].y
		Return index;
	End Method  

End Type


Type vec2
  Field x:Int,y:Int

  Function Create:vec2(x:Int,y:Int)
    Local tmp:vec2 = New vec2
    tmp.x = x; tmp.y = y;
    Return tmp
  End Function

End Type